# -*- coding: utf-8 -*-
'''
Created on May 20, 2012

@author: tel
'''
import sys, os, time, numpy, site, re
site.addsitedir('../')
from collections import namedtuple

import makedb.levels as levels
from makedb.levels_error import ConformerEmptyError
from histo import SixPanelHisto, OnePanelHisto
from scatter import Scatter, MultiScatter
from hexbin import Hexbin, MultiHexbin, MultiHexbinS
from bezier import Bezier
from sqlalchemy import Table, Column, Integer, ForeignKey, create_engine, event, MetaData
from sqlalchemy.orm import relationship, backref, aliased, sessionmaker, joinedload
from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy.pool import NullPool

def _fk_pragma_on_connect(dbapi_con, con_record):
    dbapi_con.execute('pragma foreign_keys=ON')

def regexp(expr, item):
    r = re.compile(expr)
    return r.match(item) is not None

class Stats(object):
    def __init__(self, name, eattrs=None):
        if eattrs==None:
            eattrs = []
        self.eattrs = list(eattrs)
        for eattr in self.eattrs:
            self.__setattr__(eattr, [])
        self.attrs = ['tooHigh', 'tooLow', 'highLow', 'highChi2', 'tooSharp', 'negDsol', 'badRes', 'goodRes', 'count']
        for attr in self.attrs:
            self.__setattr__(attr, 0)
        self.name = name
        self._badRes = 0
        self._goodRes = 0
        
    def Gate(self, rlist):
        self.count += 1
        hcount = 0
        lcount = 0
        hchicount = 0
        sharpcount = 0
        negdsolcount = 0
        for res in rlist:
            if res.orh==1:
                hcount += 1
            elif res.orl==1:
                lcount += 1
            elif res.toosharp==1:
                sharpcount += 1
            elif getattr(res, '1000*chi2') > 50:
                hchicount += 1
        if hcount+lcount+hchicount+sharpcount >0:
            if hcount > 0 and lcount > 0:
                self.highLow += 1
            elif hcount > 0:
                self.tooHigh += 1
            elif lcount > 0:
                self.tooLow += 1
            if hchicount > 0:
                self.highChi2 += 1
            if sharpcount > 0:
                self.tooSharp += 1
            if negdsolcount > 0:
                self.negDsol += 1
            self.badRes += 1
            return False
        else:
            self.goodRes += 1
            return True
    
    def Avg(self, attr):
        return numpy.average(self.__getattribute__(attr))
        
    def Std(self, attr):
        return numpy.std(self.__getattribute__(attr))
    
    def __str__(self):
        ps = ['%s\n' % self.name]
        for attr in self.attrs:
            ps.append("%s:\t%s\n" % (attr, self.__getattribute__(attr)))
        for eattr in self.eattrs:
            ps.append("average %s:\t%.3f\t±\t%.3f\n" % (eattr, self.Avg(eattr), self.Std(eattr)))
        return ''.join(ps)

def LoadSession(db):
    Session = sessionmaker(bind=db)
    session = Session()
    return session
        
class EpsPaperGraphs(object):
    
    class Info(object):
        dsols = ('all_dsol','low_dsol','mid_dsol','high_dsol')
        ress = ('ALL','ASP','GLU','ARG','LYS','HIS','TYR','SNASE')
        ResD = namedtuple('ResD', ('ALL','ASP','GLU','ARG','LYS','HIS','TYR','SNASE'))
        
        def __init__(self, **kwargs):
            D = namedtuple('D', ['stats']+kwargs.keys())
            for dsol in self.__class__.dsols:
                resargs = []
                for arg in self.__class__.ResD._fields:
                    resargs.append(None)
                self.__setattr__(dsol, self.__class__.ResD(*resargs))
                for res in self.__class__.ress:
                    dargs = [Stats('%s_%s' % (dsol, res))]
                    for val in kwargs.itervalues():
                        valargs = []
                        for arg in val._fields:
                            valargs.append([])
                        dargs.append(val(*valargs))
                    replacedict = {res:D(*dargs)}    
                    self.__setattr__(dsol, self.__getattribute__(dsol)._replace(**replacedict))
        
        def GenPlotInput(self, dsol, dname, x, y):
            plotinput = []
            for lis in (['ASP','GLU','ARG','LYS'], ['HIS'], ['TYR'],['SNASE']):
                for i in range(2):
                    plotinput.append([])
                for res in lis:
                    d = self.__getattribute__(dsol).__getattribute__(res).__getattribute__(dname)
                    plotinput[-2]+=d.__getattribute__(x)
                    plotinput[-1]+=d.__getattribute__(y)
            return plotinput
                        
    HistoD = namedtuple('HistoD', ['ld_dpk', 'ld_ddpk', 'md_dpk', 'md_ddpk', 'hd_dpk', 'hd_ddpk', 'mut_dpk', 'mut_ddpk'])
    ModelD = namedtuple('ModelD', ['ld_ddpk','ld_mddpk','md_ddpk','md_mddpk','hd_ddpk','hd_mddpk'])
    SizeD = namedtuple('SizeD', ['size','deepcount'])
    RMSDD = namedtuple('RMSDD', ['eps', 'wt_rmsd', 'mut_rmsd'])
    IndivD = namedtuple('IndivD', ['eps', 'res_data', 'mut_rmsd'])
    StatD = namedtuple('StatD', ['ld_stat', 'md_stat', 'hd_stat', 'ad_stat'])
    DsolFourTwentyD = namedtuple('DsolFourTwentyD', ['ld_fdsol','ld_tdsol','md_fdsol','md_tdsol','hd_fdsol','hd_tdsol'])
    EleFourTwentyD = namedtuple('EleFourTwentyD', ['ld_fele','ld_tele','md_fele','md_tele','hd_fele','hd_tele'])
    RatioFourTwentyD = namedtuple('RatioFourTwentyD', ['ld_fratio','ld_tratio','md_fratio','md_tratio','hd_fratio','hd_tratio'])
    QFdpkD = namedtuple('QFdpkD', ['ld_qdpk', 'ld_fdpk', 'md_qdpk', 'md_fdpk', 'hd_qdpk', 'hd_fdpk'])
    QFddpkD = namedtuple('QFddpkD', ['ld_qddpk', 'ld_fddpk', 'md_qddpk', 'md_fddpk', 'hd_qddpk', 'hd_fddpk'])
    SetsDsolEleD = namedtuple('DsolEleD', ['dsolderk', 'elederk', 'dsolhistyr', 'elehistyr', 'dsolsnase', 'elesnase'])
    houk_res = (('3RN3','ASP',14,1.8),('4LZT','ASP',87,2.1),('1PPF','ASP',27,2.2),('1XNB','ASP',11,2.5),('1BEO','ASP',21,2.5),('4LZT','ASP',18,2.7),('1XNB','ASP',106,2.7),('1PGA','ASP',22,2.9),('3RN3','ASP',121,3.1),('1A2P','ASP',75,3.1),('2RN2','ASP',94,3.2),('1PGA','ASP',47,3.4),('3RN3','ASP',53,3.7),('4LZT','ASP',52,3.7),('1PGA','ASP',36,3.8),('2TRX','ASP',20,3.8),('1DE3','ASP',59,4.1),('1DE3','ASP',57,4.3),('2RN2','ASP',10,6.1),('2TRX','ASP',26,8.1),('3RN3','GLU',2,2.6),('4LZT','GLU',7,2.9),('1PPF','GLU',19,3.2),('2RN2','GLU',57,3.2),('1A2P','GLU',60,3.2),('3RN3','GLU',111,3.5),('2RN2','GLU',129,3.6),('2RN2','GLU',61,3.9),('3RN3','GLU',9,4),('1BCA','GLU',26,4.1),('1PPF','GLU',10,4.1),('2RN2','GLU',119,4.1),('1PGA','GLU',27,4.5),('1PPF','GLU',43,4.8),('1DE3','GLU',96,5.1),('1ANS','GLU',20,5.4),('1RGA','GLU',28,5.9),('4LZT','GLU',35,6.2),('1HNG','GLU',41,6.7),('1XNB','GLU',172,6.7),('2SNM','LYS',66,6.4),('1L54','LYS',102,6.6),('1MUT','LYS',39,8.4),('1NFN','LYS',146,9.2),('1FEZ','LYS',53,9.3),('1GS9','LYS',146,9.4),('1LE2','LYS',143,9.4),('1NFN','LYS',143,9.5),('1NZP','LYS',312,9.5),('1GS9','LYS',143,9.9),('1LE2','LYS',146,9.9),('1PPF','LYS',34,10.1),('4LZT','LYS',33,10.4),('2BCA','LYS',41,10.8),('4LZT','LYS',96,10.8),('1PGA','LYS',28,10.9),('2BCA','LYS',16,10.9),('1PPF','LYS',55,11.1),('2BCA','LYS',7,11.2),('2BCA','LYS',55,12.4),('3EBX','HIS',6,2.8),('3SSI','HIS',43,3.2),('1STN','HIS',121,5.3),('4LZT','HIS',15,5.4),('1ERT','HIS',43,5.5),('1DE3','HIS',137,5.8),('3RN3','HIS',48,6.1),('3RN3','HIS',119,6.5),('3RN3','HIS',12,6),('1DE3','HIS',104,6.5),('1DE3','HIS',36,6.8),('2RN2','HIS',62,7),('2RN2','HIS',124,7.1),('1DE3','HIS',50,7.7),('1RGA','HIS',92,7.8),('1RGA','HIS',40,7.9),('2RN2','HIS',127,7.9),('1DG9','HIS',66,8.3),('2LZM','HIS',31,9.1),('1DG9','HIS',72,9.2))   
    in_silico_res = [ [ 'V023D', 'ASP', '23', '6.8' ], [ 'L025D', 'ASP', '25', '6.8' ], [ 'F034D', 'ASP', '34', '7.8' ], [ 'L036D', 'ASP', '36', '7.9' ], [ 'L037D', 'ASP', '37', '<4.0' ], [ 'L038D', 'ASP', '38', '6.8' ], [ 'V039D', 'ASP', '39', '8.1' ], [ 'T041D', 'ASP', '41', '<4.0' ], [ 'A058D', 'ASP', '58', '6.8' ], [ 'T062D', 'ASP', '62', '8.7' ], [ 'V066D', 'ASP', '66', '8.7' ], [ 'I072D', 'ASP', '72', '7.6' ], [ 'V074D', 'ASP', '74', '8.3' ], [ 'A090D', 'ASP', '90', '7.5' ], [ 'Y091D', 'ASP', '91', '7.2' ], [ 'I092D', 'ASP', '92', '8.1' ], [ 'V099D', 'ASP', '99', '8.5' ], [ 'N100D', 'ASP', '100', '6.9' ], [ 'L103D', 'ASP', '103', '8.7' ], [ 'V104D', 'ASP', '104', '9.7' ], [ 'A109D', 'ASP', '109', '7.5' ], [ 'N118D', 'ASP', '118', '7.0' ], [ 'L125D', 'ASP', '125', '7.6' ], [ 'A132D', 'ASP', '132', '7.0' ], [ 'V023E', 'GLU', '23', '7.1' ], [ 'L025E', 'GLU', '25', '7.5' ], [ 'F034E', 'GLU', '34', '7.3' ], [ 'L036E', 'GLU', '36', '8.7' ], [ 'L037E', 'GLU', '37', '5.2' ], [ 'L038E', 'GLU', '38', '7.0' ], [ 'V039E', 'GLU', '39', '8.2' ], [ 'T041E', 'GLU', '41', '6.5' ], [ 'A058E', 'GLU', '58', '7.7' ], [ 'T062E', 'GLU', '62', '7.7' ], [ 'V066E', 'GLU', '66', '8.5' ], [ 'I072E', 'GLU', '72', '7.3' ], [ 'V074E', 'GLU', '74', '7.8' ], [ 'A090E', 'GLU', '90', '6.4' ], [ 'Y091E', 'GLU', '91', '7.3' ], [ 'I092E', 'GLU', '92', '9.0' ], [ 'V099E', 'GLU', '99', '8.4' ], [ 'N100E', 'GLU', '100', '7.6' ], [ 'L103E', 'GLU', '103', '8.9' ], [ 'V104E', 'GLU', '104', '9.4' ], [ 'A109E', 'GLU', '109', '7.9' ], [ 'N118E', 'GLU', '118', '<4.5' ], [ 'L125E', 'GLU', '125', '9.1' ], [ 'A132E', 'GLU', '132', '7.0' ], [ 'V023K', 'LYS', '23', '7.3' ], [ 'L025K', 'LYS', '25', '6.3' ], [ 'F034K', 'LYS', '34', '7.1' ], [ 'L036K', 'LYS', '36', '7.2' ], [ 'L037K', 'LYS', '37', '>10.4' ], [ 'L038K', 'LYS', '38', '>10.4' ], [ 'V039K', 'LYS', '39', '9.0' ], [ 'T041K', 'LYS', '41', '9.3' ], [ 'A058K', 'LYS', '58', '10.4' ], [ 'T062K', 'LYS', '62', '8.1' ], [ 'V066K', 'LYS', '66', '5.6' ], [ 'I072K', 'LYS', '72', '8.6' ], [ 'V074K', 'LYS', '74', '7.4' ], [ 'A090K', 'LYS', '90', '8.6' ], [ 'Y091K', 'LYS', '91', '5.3' ], [ 'I092K', 'LYS', '92', '5.3' ], [ 'V099K', 'LYS', '99', '6.5' ], [ 'N100K', 'LYS', '100', '8.6' ], [ 'L103K', 'LYS', '103', '8.2' ], [ 'V104K', 'LYS', '104', '7.7' ], [ 'A109K', 'LYS', '109', '9.2' ], [ 'N118K', 'LYS', '118', '>10.4' ], [ 'L125K', 'LYS', '125', '6.2' ], [ 'A132K', 'LYS', '132', '>10.4' ], [ 'V023R', 'ARG', '23', '>10.4' ], [ 'L025R', 'ARG', '25', '>10.4' ], [ 'F034R', 'ARG', '34', '>10.4' ], [ 'L036R', 'ARG', '36', '>10.4' ], [ 'L037R', 'ARG', '37', '>10.4' ], [ 'L038R', 'ARG', '38', '>10.4' ], [ 'V039R', 'ARG', '39', '>10.4' ], [ 'T041R', 'ARG', '41', '>10.4' ], [ 'A058R', 'ARG', '58', '>10.4' ], [ 'T062R', 'ARG', '62', '>10.4' ], [ 'V066R', 'ARG', '66', '>10.4' ], [ 'I072R', 'ARG', '72', '>10.4' ], [ 'V074R', 'ARG', '74', '>10.4' ], [ 'A090R', 'ARG', '90', '>10.4' ], [ 'Y091R', 'ARG', '91', '>10.4' ], [ 'I092R', 'ARG', '92', '>10.4' ], [ 'V099R', 'ARG', '99', '>10.4' ], [ 'N100R', 'ARG', '100', '>10.4' ], [ 'L103R', 'ARG', '103', '>10.4' ], [ 'V104R', 'ARG', '104', '>10.4' ], [ 'A109R', 'ARG', '109', '>10.4' ], [ 'N118R', 'ARG', '118', '>10.4' ], [ 'L125R', 'ARG', '125', '>10.4' ], [ 'A132R', 'ARG', '132', '>10.4' ] ]
    
    DsolEleD = namedtuple('DsolEleD', ['dsol', 'ele', 'asite_dist'])
    RatioD = namedtuple('RatioD', ['ratio', 'ddpk'])
    DpkDDpkD = namedtuple('DpkDDpkD', ['dpk', 'ddpk'])
    
    def __init__(self, pdbsession, shsession, insilicosession, csasession=None, consurfsession=None, outdir='/home/tel/benchmark/links'):
        temp = []
        for res in self.__class__.in_silico_res:
            if int(res[2]) < 50:
                res[2] = int(res[2]) - 6
            else:
                res[2] = int(res[2]) - 12
            #res[3] = float(res[3])
            temp.append(res)
        self.in_silico_res = temp
        self.dsoldict = {'all_dsol':(-999,999),'low_dsol':(-999,3),'mid_dsol':(3,6),'high_dsol':(6,999)}
        
        self.pdbsession = pdbsession
        self.shsession = shsession
        self.insilicosession = insilicosession
        self.csasession = csasession
        self.consurfsession = consurfsession
        self.outdir = outdir
        self.histod = self.__class__.HistoD([], [], [], [], [], [],
                                            [10.4,10.4,10.4,10.4,10.4,10.4,10.4,10.4,10.4,10.4,9.25,9.25,9.043,10.4,10.4,9.25,9.25,9.25,9.25,9.25,9.25,10.4,5.274,9.25,10.4,8.567,9.25,10.4,9.25,8.602,9.25,9.25,10.4,10.4,8.572,10.4,10.189,9.25,10.4,7.337,9.25,7.234,9.25,6.838,7.899,7.271,6.568,8.861,9.25,10.4,9.25,5.729,9.25,7.647,7.274,7.282,6.627,7.603,5.959,9.25,9.25,1.18,9.25,9.198,7.87],
                                            [-8.621,-8.589,-9.087,-9.717,-9.804,-9.645,-8.461,-10.468,-9.336,-8.953,-8.502,-8.541,-8.529,-9.375,-8.628,-7.72,-7.738,-8.143,-8.067,-7.774,-8.356,-8.714,-6.201,-7.433,-8.303,-8.202,-8.767,-9.755,-8.652,-7.954,-8.381,-8.294,-11.012,-9.544,-8.115,-8.537,-11.551,-8.532,-12.594,-7.095,-7.545,-7.416,-8.231,-6.486,-7.255,-6.653,-6.043,-8.54,-8.211,-8.672,-8.246,-5.19,-8.948,-7.002,-6.486,-6.584,-6.019,-7.275,-5.741,-8.188,-9.437,-3.061,-6.741,-7.432,-6.6])
        self.ratiod = self.__class__.RatioD([],[])
        self.modeld = self.__class__.ModelD([],[],[],[],[],[])
        self.sized = self.__class__.SizeD([],[])
        self.dsoleled = self.__class__.DsolEleD([],[],[])
        self.qfdpkd= self.__class__.QFdpkD([],[],[],[],[],[])
        self.qfddpkd= self.__class__.QFddpkD([],[],[],[],[],[])
        self.dsolfourtwentyd= self.__class__.DsolFourTwentyD([],[],[],[],[],[])
        self.elefourtwentyd= self.__class__.EleFourTwentyD([],[],[],[],[],[])
        self.ratiofourtwentyd= self.__class__.RatioFourTwentyD([],[],[],[],[],[])
        self.rmsdd = self.__class__.RMSDD([2,4,6,8,10,12,20],[],[7.09,6.36,4.34,2.18,1.68,1.68,2.39])   #list of mutant RMSDs with arginine included [7.80722355,5.978183016,3.860512303,1.914451793,1.455808442,1.451851989,2.065739655]
        self.indivd = self.__class__.IndivD([2,4,6,8,10,12,20],[],[7.09,6.36,4.34,2.18,1.68,1.68,2.39])
        self.statd = self.__class__.StatD(Stats('Low Dsol Stats'), Stats('Mid Dsol Stats'), Stats('High Dsol Stats'), Stats('All Dsol Stats'))
        
        self.infokwargs = {}
        for key in ['DsolEleD', 'RatioD', 'DpkDDpkD']:
            self.infokwargs[key.lower()] = self.__getattribute__(key)
        self.info = self.__class__.Info(**self.infokwargs)
            
    def GoSHAggregate(self):
        total = time.clock()
        now = time.clock()
        count = 0
        
        for eps in (2,4,6,8,10,12,20):
            tmp = []
            for pdb, res, resn, pk in self.__class__.houk_res:
                shres = self.shsession.query(levels.ResData).join(levels.RunData, levels.ProtData).\
                             filter(levels.ProtData.pdb==pdb).filter(levels.RunData.EPSILON_PROT==eps).\
                             filter(levels.ResData.res==res).filter(levels.ResData.resn==resn).\
                             filter(levels.ResData.seqadv=='').first()
                if shres!=None:
                    if shres.res!='HIS':
                        tmp.append((pk - shres.__getattribute__('pKa/Em'))**2)
                        #print shres, ('\t%s\t%s\t%s' % (shres.e_r, shres.dsol, shres.dpk))
            self.rmsdd.wt_rmsd.append(numpy.sqrt(numpy.mean(tmp)))
            print '%d\tdielectric constant %d took: %.3f' % (count, eps, time.clock() - now)
            now = time.clock()
            count+=1
        Bezier(*self.rmsdd, outdir=self.outdir, outfile='rmsd_bezier')
        print 'the whole thing took: %.3f' % (time.clock() - total)
    
    def GoSHIndividual(self):
        total = time.clock()
        now = time.clock()
        count = 0
        epss = [2,4,6,8,10,12,20]
        
        for pdb, res, resn, pk in self.__class__.houk_res:
            tmp = []
            for eps in epss:
                shres = self.shsession.query(levels.ResData).join(levels.RunData, levels.ProtData).\
                             filter(levels.ProtData.pdb==pdb).filter(levels.RunData.EPSILON_PROT==eps).\
                             filter(levels.ResData.res==res).filter(levels.ResData.resn==resn).\
                             filter(levels.ResData.seqadv=='').first()
                if shres!=None:
                    if shres.res!='HIS':
                        tmp.append(numpy.absolute(pk - shres.__getattribute__('pKa/Em')))
                        #print shres, ('\t%s\t%s\t%s' % (shres.e_r, shres.dsol, shres.dpk))
            self.indivd.res_data.append(tmp)
            print '%d\tdielectric constant %d took: %.3f' % (count, eps, time.clock() - now)
            now = time.clock()
            count+=1
        Bezier(*self.indivd, outdir=self.outdir, outfile='rmsd_indiv_bezier')
        print 'the whole thing took: %.3f' % (time.clock() - total)
    
    def GoSets(self):
        total = time.clock()
        now = time.clock()
        count = 0
        for prot in self.pdbsession.query(levels.ProtData).yield_per(10):
            rundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
            self.KerchunkSets(rundict, snase=False)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        for prot in self.insilicosession.query(levels.ProtData).yield_per(10):
            rundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
            self.KerchunkSets(rundict, snase=True)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        
#        x1 = self.info.all_dsol.ASP.dsoleled.dsol + self.info.all_dsol.GLU.dsoleled.dsol + self.info.all_dsol.ARG.dsoleled.dsol + self.info.all_dsol.LYS.dsoleled.dsol
#        y1 = self.info.all_dsol.ASP.dsoleled.ele + self.info.all_dsol.GLU.dsoleled.ele + self.info.all_dsol.ARG.dsoleled.ele + self.info.all_dsol.LYS.dsoleled.ele
#        x2 = self.info.all_dsol.HIS.dsoleled.dsol + self.info.all_dsol.TYR.dsoleled.dsol
#        y2 = self.info.all_dsol.HIS.dsoleled.ele + self.info.all_dsol.TYR.dsoleled.ele
#        x3 = self.info.all_dsol.SNASE.dsoleled.dsol
#        y3 = self.info.all_dsol.SNASE.dsoleled.ele
        
        print 'graphing now...'
        MultiHexbin(*self.info.GenPlotInput(dsol='all_dsol', dname='dsoleled', x='dsol', y='ele'), outdir=self.outdir, outfile='dsol_ele-alldsol-derk_h_y_snase-multihexbin')
        MultiHexbin(*self.info.GenPlotInput(dsol='high_dsol', dname='ratiod', x='ratio', y='ddpk'), outdir=self.outdir, outfile='ratio_ddpk-highdsol-derk_h_y_snase-multihexbin')
        derk_ld = self.info.GenPlotInput(dsol='low_dsol', dname='dpkddpkd', x='dpk', y='ddpk')[0:2]
        derk_md = self.info.GenPlotInput(dsol='mid_dsol', dname='dpkddpkd', x='dpk', y='ddpk')[0:2]
        derk_hd = self.info.GenPlotInput(dsol='high_dsol', dname='dpkddpkd', x='dpk', y='ddpk')[0:2]
        his_ld = self.info.GenPlotInput(dsol='low_dsol', dname='dpkddpkd', x='dpk', y='ddpk')[2:4]
        his_md = self.info.GenPlotInput(dsol='mid_dsol', dname='dpkddpkd', x='dpk', y='ddpk')[2:4]
        his_hd = self.info.GenPlotInput(dsol='high_dsol', dname='dpkddpkd', x='dpk', y='ddpk')[2:4]
        snase_hd = self.info.GenPlotInput(dsol='high_dsol', dname='dpkddpkd', x='dpk', y='ddpk')[6:8]
        SixPanelHisto(*(derk_ld+derk_md+derk_hd+snase_hd), outdir=self.outdir, outfile='dpka_ddpka-lowdsol_middsol_highdsol-derk-histos')
        SixPanelHisto(*(his_ld+his_md+his_hd+snase_hd), outdir=self.outdir, outfile='dpka_ddpka-lowdsol_middsol_highdsol-his-histos')
        print 'the whole thing took: %.3f' % (time.clock() - total)
            
    def GoDun(self):
        total = time.clock()
        now = time.clock()
        count = 0
        for prot in self.pdbsession.query(levels.ProtData).yield_per(10):
            rundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
            self.Kerchunk(rundict)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        
        print 'graphing now...'
        MultiScatter(*self.dsolfourtwentyd, outdir=self.outdir, outfile='e4_dsol_e20_dsol_multiscatter')
        MultiScatter(*self.elefourtwentyd, outdir=self.outdir, outfile='e4_ele_e20_ele_multiscatter')
        MultiScatter(*self.ratiofourtwentyd, outdir=self.outdir, outfile='e4_ratio_e20_ratio_multiscatter')
        Scatter(self.modeld[4], self.modeld[5], outdir=self.outdir, outfile='ddpk_mddpk_scatter')
        MultiHexbin(*self.dsolfourtwentyd, outdir=self.outdir, outfile='e4-dsol_e20-dsol_multihexbin')
        MultiHexbin(*self.elefourtwentyd, outdir=self.outdir, outfile='e4-ele_e20-ele_multihexbin')
        MultiHexbin(*self.ratiofourtwentyd, outdir=self.outdir, outfile='e4-ratio_e20-ratio_multihexbin')
        Scatter(*self.ratiod, outdir=self.outdir, outfile='ratio_ddpk_scatter')
        Scatter(self.modeld[4], self.modeld[5], outdir=self.outdir, outfile='ddpk_mddpk_scatter')
        Hexbin(*self.dsoleled, outdir=self.outdir, outfile='dsol_ele_hexbin')
        Hexbin(*self.ratiod, outdir=self.outdir, outfile='ratio_ddpk_hexbin')
        SixPanelHisto(*self.histod, outdir=self.outdir, outfile='dpka_ddpka_histos')
        print 'the whole thing took: %.3f' % (time.clock() - total)
        
    def GoDunFQ(self):
        total = time.clock()
        now = time.clock()
        count = 0
        for prot in self.pdbsession.query(levels.ProtData).yield_per(10):
            rundict = {}
            fullrundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
                else:
                    fullrundict[run.EPSILON_PROT] = run
            self.KerchunkFQ(rundict, fullrundict)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        
        print 'graphing now...'
        MultiScatter(*self.qfdpkd, outdir=self.outdir, outfile='quick_dpk_full_dpk_multiscatter')
        MultiScatter(*self.qfddpkd, outdir=self.outdir, outfile='quick_ddpk_full_ddpk_multiscatter')
        print 'done scatter'
        MultiHexbin(*self.qfdpkd, outdir=self.outdir, outfile='quick_dpk_full_dpk_multihexbin')
        MultiHexbin(*self.qfddpkd, outdir=self.outdir, outfile='quick_ddpk_full_ddpk_multihexbin')
        print 'the whole thing took: %.3f' % (time.clock() - total)

    def GoDunCSA(self):
        total = time.clock()
        now = time.clock()
        count = 0
        self.csacount = 0
        for prot in self.pdbsession.query(levels.ProtData).yield_per(10):
            rundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
            self.KerchunkCSA(rundict)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        
        print 'graphing now...'
        Scatter(*self.dsoleled, outdir=self.outdir, outfile='dsol_ele_csa_scatter')
        print 'the whole thing took: %.3f' % (time.clock() - total)
    
    def GoProtStat(self):
        total = time.clock()
        now = time.clock()
        count = 0
        self.csacount = 0
        for prot in self.pdbsession.query(levels.ProtData).yield_per(10):
            rundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
            if 4 in rundict and 20 in rundict:
                self.sized.size.append(len(rundict[4].Ress)-2)    #the -2 takes care of the added terminals: NTR and CTR
                deepcount = 0
                for res4 in rundict[4].Ress:
                    if res4.dsol > 6:
                        deepcount+=1
                self.sized.deepcount.append(deepcount)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        
        print 'graphing now...'
        Hexbin(*self.sized, outdir=self.outdir, outfile='size_deepcount_hexbin')
        print '%s\t%s\t%s\t%s' % ('residues >', 'residues <=', 'total proteins', 'highly buried res')
        for l,r in ((0,200),(200,400),(400,600),(600,800),(800,1000),(1000,max(self.sized.size))):
            tempdcount = []
            total = 0
            for size, deepcount in zip(*self.sized):
                if l < size <= r:
                    total += 1
                    tempdcount.append(deepcount)
            if tempdcount!=[]:
                print '%d\t%d\t%d\t%3.3f\t+-\t%3.3f' % (l,r,total,numpy.average(tempdcount),numpy.std(tempdcount))
            else:
                print '%d\t%d\t%d\t%s\t+-\t%s' % (l,r,total,'none','none')
        print 'the whole thing took: %.3f' % (time.clock() - total)
    
    def Kerchunk(self, rundict):
        if 4 in rundict and 20 in rundict:
            for res4, res20 in zip(rundict[4].Ress, rundict[20].Ress):
                if res4==res20 and res4.res in ['ASP','GLU','ARG','LYS'] and res4.dpk!=None and res4.seqadv=='':
                    if res4.dsol <= 3 and self.statd.ld_stat.Gate([res4, res20]):
                        dsol4 = 0.001 if res4.dsol <= 0 else res4.dsol
                        dsol20 = 0.001 if res20.dsol <= 0 else res20.dsol
                        self.histod.ld_dpk.append(res4.dpk)
                        self.histod.ld_ddpk.append(res20.dpk - res4.dpk)
                        self.dsolfourtwentyd.ld_fdsol.append(dsol4)
                        self.dsolfourtwentyd.ld_tdsol.append(dsol20)
                        self.elefourtwentyd.ld_fele.append(res4.e_r + res4.e_b)
                        self.elefourtwentyd.ld_tele.append(res20.e_r + res20.e_b)
                        if res4.dsol>.5:
                            self.modeld.hd_mddpk.append((16*((res4.e_r + res4.e_b)*(76)+res4.dsol*80))/(20*(-76)))
                            self.ratiofourtwentyd.ld_fratio.append((res4.e_r + res4.e_b)/(dsol4))
                            self.ratiofourtwentyd.ld_tratio.append((res20.e_r + res20.e_b)/(dsol20))
                    elif 3 < res4.dsol <= 6 and self.statd.md_stat.Gate([res4, res20]):
                        self.histod.md_dpk.append(res4.dpk)
                        self.histod.md_ddpk.append(res20.dpk - res4.dpk)
                        self.dsolfourtwentyd.md_fdsol.append(res4.dsol)
                        self.dsolfourtwentyd.md_tdsol.append(res20.dsol)
                        self.elefourtwentyd.md_fele.append(res4.e_r + res4.e_b)
                        self.elefourtwentyd.md_tele.append(res20.e_r + res20.e_b)
                        self.modeld.md_mddpk.append((16*((res4.e_r + res4.e_b)*(76)+res4.dsol*80))/(20*(-76)))
                        if res4.dsol>.5:
                            self.ratiofourtwentyd.md_fratio.append((res4.e_r + res4.e_b)/(res4.dsol))
                            self.ratiofourtwentyd.md_tratio.append((res20.e_r + res20.e_b)/(res20.dsol))
                    elif 6 < res4.dsol and self.statd.hd_stat.Gate([res4, res20]):
                        self.histod.hd_dpk.append(res4.dpk)
                        for l in (self.histod.hd_ddpk,self.ratiod.ddpk,self.modeld.hd_ddpk):
                            l.append(res20.dpk - res4.dpk)
                        self.ratiod.ratio.append(float(res4.e_r + res4.e_b)/res4.dsol)
                        self.modeld.hd_mddpk.append((16*((res4.e_r + res4.e_b)*(76)+res4.dsol*80))/(20*(-76)))
                        self.dsolfourtwentyd.hd_fdsol.append(res4.dsol)
                        self.dsolfourtwentyd.hd_tdsol.append(res20.dsol)
                        self.elefourtwentyd.hd_fele.append(res4.e_r + res4.e_b)
                        self.elefourtwentyd.hd_tele.append(res20.e_r + res20.e_b)
                        if res4.dsol>.5:
                            self.ratiofourtwentyd.hd_fratio.append((res4.e_r + res4.e_b)/(res4.dsol))
                            self.ratiofourtwentyd.hd_tratio.append((res20.e_r + res20.e_b)/(res20.dsol))
                    if self.statd.ad_stat.Gate([res4, res20]):
                        self.dsoleled.ele.append(res4.e_r + res4.e_b)
                        self.dsoleled.dsol.append(res4.dsol)
                else:
                    pass
                    #print "(%s %s %s) (%s %s %s)" % (res4.res, res4.chainid, res4.resn, res20.res, res20.chainid, res20.resn)
    def KerchunkSets(self, rundict, snase):
        if 4 in rundict and 20 in rundict:
            for res4, res20 in zip(rundict[4].Ress, rundict[20].Ress):
                if snase==True:
                    if (res4.Prot.pdb, res4.res, res4.resn) in zip(*zip(*self.in_silico_res)[0:3]):
                        pass
                    else:
                        continue
                if res4.dpk!=None and res4.seqadv=='':
                    if snase==True:
                        infolist = ['SNASE']
                    elif res4.res in ['ASP','GLU','ARG','LYS']:  
                        infolist = [res4.res,'ALL']
                    elif res4.res in ['HIS','TYR']:
                        infolist = [res4.res]
                    else:
                        infolist = []
                else:
                    continue
                for dsol, dsolrange in self.dsoldict.iteritems():
                    if (dsolrange[0] < res4.dsol <= dsolrange[1]):
                        for r in infolist:
                            info = self.info.__getattribute__(dsol).__getattribute__(r)
                            if info.stats.Gate([res4,res20]):
                                if res4.dsol<=0:
                                    dsol4 = .001
                                else:
                                    dsol4 = res4.dsol
                                info.dsoleled.dsol.append(dsol4)
                                info.dsoleled.ele.append(res4.e_b + res4.e_r)
                                info.dpkddpkd.dpk.append(res4.dpk)
                                info.dpkddpkd.ddpk.append(res20.dpk - res4.dpk)
                                if dsol4 > .5:
                                    info.ratiod.ratio.append((res4.e_b + res4.e_r)/dsol4)
                                    info.ratiod.ddpk.append(res20.dpk - res4.dpk)
    
    def KerchunkFQ(self, rundict, fullrundict):
        '''
        kerchunker for full vs quick analysis
        '''
        if 4 in rundict and 20 in rundict and 4 in fullrundict and 20 in fullrundict:
            for qres4, qres20, fres4, fres20 in zip(rundict[4].Ress, rundict[20].Ress, fullrundict[4].Ress, fullrundict[20].Ress):
                if qres4==qres20==fres4==fres20 and qres4.res in ['ASP','GLU','HIS','ARG','LYS'] and qres4.dpk!=None and qres4.seqadv=='':
                    if qres4.dsol <= 3 and self.statd.ld_stat.Gate([qres4, qres20, fres4, fres20]):
                        self.qfdpkd.ld_qdpk.append(qres4.dpk)
                        self.qfdpkd.ld_fdpk.append(fres4.dpk)
                        self.qfddpkd.ld_qddpk.append(qres20.dpk - qres4.dpk)
                        self.qfddpkd.ld_fddpk.append(fres20.dpk - fres4.dpk)
                    elif 3 < qres4.dsol <= 6 and self.statd.md_stat.Gate([qres4, qres20, fres4, fres20]):
                        self.qfdpkd.md_qdpk.append(qres4.dpk)
                        self.qfdpkd.md_fdpk.append(fres4.dpk)
                        self.qfddpkd.md_qddpk.append(qres20.dpk - qres4.dpk)
                        self.qfddpkd.md_fddpk.append(fres20.dpk - fres4.dpk)
                    elif 6 < qres4.dsol and self.statd.hd_stat.Gate([qres4, qres20, fres4, fres20]):
                        self.qfdpkd.hd_qdpk.append(qres4.dpk)
                        self.qfdpkd.hd_fdpk.append(fres4.dpk)
                        self.qfddpkd.hd_qddpk.append(qres20.dpk - qres4.dpk)
                        self.qfddpkd.hd_fddpk.append(fres20.dpk - fres4.dpk)
                else:
                    pass
                    #print "(%s %s %s) (%s %s %s)" % (res4.res, res4.chainid, res4.resn, res20.res, res20.chainid, res20.resn)
    
    def KerchunkCSA(self, rundict):
        if 4 in rundict and 20 in rundict:
            actress = []
            try:
                csas = self.csasession.query(levels.CSAData).filter(levels.CSAData.pdb==rundict[4].Ress[0].Prot.pdb.lower()).all()
                for csa in csas:
                    if len(csa.res) < 3:
                        #needs to be csares instead of csa.res to avoid updating the database
                        csares = '_'+csa.res
                    elif len(csa.res) < 1:
                        continue
                    else:
                        csares = csa.res
                    actres = self.pdbsession.query(levels.ResData).join(levels.ProtData).\
                    filter(levels.ProtData.pdb==csa.pdb.upper()).filter(levels.ResData.res==csares).\
                    filter(levels.ResData.resn==csa.resn).filter(levels.ResData.chainid==csa.chainid).first()
                    if actres==None and csares=='CYS':
                        actres = self.pdbsession.query(levels.ResData).join(levels.ProtData).\
                        filter(levels.ProtData.pdb==csa.pdb.upper()).filter(levels.ResData.res=='CYD').\
                        filter(levels.ResData.resn==csa.resn).filter(levels.ResData.chainid==csa.chainid).first()
                    if actres!=None:
                        actress.append(actres)
            except IndexError:
                pass
            if actress!=[]:
                self.csacount += 1
            for res4, res20 in zip(rundict[4].Ress, rundict[20].Ress):
                if res4==res20 and res4.res in ['ASP','GLU','HIS','ARG','LYS'] and res4.dpk!=None and res4.seqadv=='':
                    if self.statd.hd_stat.Gate([res4, res20]):
                        self.dsoleled.ele.append(res4.e_r + res4.e_b)
                        self.dsoleled.dsol.append(res4.dsol)
                        mini = 999
                        for actres in actress:
                            try:
                                temp_min = res4.MinDist(actres)
                            except (ConformerEmptyError, AttributeError):
                                temp_min = 999
                            if temp_min < min:
                                mini = temp_min
                        self.dsoleled.asite_dist.append(mini)
                    else:
                        self.dsoleled.asite_dist.append(999)
                else:
                    pass
                    #print "(%s %s %s) (%s %s %s)" % (res4.res, res4.chainid, res4.resn, res20.res, res20.chainid, res20.resn)
    
        
    def Process(self):
        for stat in self.statd:
            print stat

class Averages(object):
    in_silico_res = [ [ 'V023D', 'ASP', '23', '6.8' ], [ 'L025D', 'ASP', '25', '6.8' ], [ 'F034D', 'ASP', '34', '7.8' ], [ 'L036D', 'ASP', '36', '7.9' ], [ 'L037D', 'ASP', '37', '<4.0' ], [ 'L038D', 'ASP', '38', '6.8' ], [ 'V039D', 'ASP', '39', '8.1' ], [ 'T041D', 'ASP', '41', '<4.0' ], [ 'A058D', 'ASP', '58', '6.8' ], [ 'T062D', 'ASP', '62', '8.7' ], [ 'V066D', 'ASP', '66', '8.7' ], [ 'I072D', 'ASP', '72', '7.6' ], [ 'V074D', 'ASP', '74', '8.3' ], [ 'A090D', 'ASP', '90', '7.5' ], [ 'Y091D', 'ASP', '91', '7.2' ], [ 'I092D', 'ASP', '92', '8.1' ], [ 'V099D', 'ASP', '99', '8.5' ], [ 'N100D', 'ASP', '100', '6.9' ], [ 'L103D', 'ASP', '103', '8.7' ], [ 'V104D', 'ASP', '104', '9.7' ], [ 'A109D', 'ASP', '109', '7.5' ], [ 'N118D', 'ASP', '118', '7.0' ], [ 'L125D', 'ASP', '125', '7.6' ], [ 'A132D', 'ASP', '132', '7.0' ], [ 'V023E', 'GLU', '23', '7.1' ], [ 'L025E', 'GLU', '25', '7.5' ], [ 'F034E', 'GLU', '34', '7.3' ], [ 'L036E', 'GLU', '36', '8.7' ], [ 'L037E', 'GLU', '37', '5.2' ], [ 'L038E', 'GLU', '38', '7.0' ], [ 'V039E', 'GLU', '39', '8.2' ], [ 'T041E', 'GLU', '41', '6.5' ], [ 'A058E', 'GLU', '58', '7.7' ], [ 'T062E', 'GLU', '62', '7.7' ], [ 'V066E', 'GLU', '66', '8.5' ], [ 'I072E', 'GLU', '72', '7.3' ], [ 'V074E', 'GLU', '74', '7.8' ], [ 'A090E', 'GLU', '90', '6.4' ], [ 'Y091E', 'GLU', '91', '7.3' ], [ 'I092E', 'GLU', '92', '9.0' ], [ 'V099E', 'GLU', '99', '8.4' ], [ 'N100E', 'GLU', '100', '7.6' ], [ 'L103E', 'GLU', '103', '8.9' ], [ 'V104E', 'GLU', '104', '9.4' ], [ 'A109E', 'GLU', '109', '7.9' ], [ 'N118E', 'GLU', '118', '<4.5' ], [ 'L125E', 'GLU', '125', '9.1' ], [ 'A132E', 'GLU', '132', '7.0' ], [ 'V023K', 'LYS', '23', '7.3' ], [ 'L025K', 'LYS', '25', '6.3' ], [ 'F034K', 'LYS', '34', '7.1' ], [ 'L036K', 'LYS', '36', '7.2' ], [ 'L037K', 'LYS', '37', '>10.4' ], [ 'L038K', 'LYS', '38', '>10.4' ], [ 'V039K', 'LYS', '39', '9.0' ], [ 'T041K', 'LYS', '41', '9.3' ], [ 'A058K', 'LYS', '58', '10.4' ], [ 'T062K', 'LYS', '62', '8.1' ], [ 'V066K', 'LYS', '66', '5.6' ], [ 'I072K', 'LYS', '72', '8.6' ], [ 'V074K', 'LYS', '74', '7.4' ], [ 'A090K', 'LYS', '90', '8.6' ], [ 'Y091K', 'LYS', '91', '5.3' ], [ 'I092K', 'LYS', '92', '5.3' ], [ 'V099K', 'LYS', '99', '6.5' ], [ 'N100K', 'LYS', '100', '8.6' ], [ 'L103K', 'LYS', '103', '8.2' ], [ 'V104K', 'LYS', '104', '7.7' ], [ 'A109K', 'LYS', '109', '9.2' ], [ 'N118K', 'LYS', '118', '>10.4' ], [ 'L125K', 'LYS', '125', '6.2' ], [ 'A132K', 'LYS', '132', '>10.4' ], [ 'V023R', 'ARG', '23', '>10.4' ], [ 'L025R', 'ARG', '25', '>10.4' ], [ 'F034R', 'ARG', '34', '>10.4' ], [ 'L036R', 'ARG', '36', '>10.4' ], [ 'L037R', 'ARG', '37', '>10.4' ], [ 'L038R', 'ARG', '38', '>10.4' ], [ 'V039R', 'ARG', '39', '>10.4' ], [ 'T041R', 'ARG', '41', '>10.4' ], [ 'A058R', 'ARG', '58', '>10.4' ], [ 'T062R', 'ARG', '62', '>10.4' ], [ 'V066R', 'ARG', '66', '>10.4' ], [ 'I072R', 'ARG', '72', '>10.4' ], [ 'V074R', 'ARG', '74', '>10.4' ], [ 'A090R', 'ARG', '90', '>10.4' ], [ 'Y091R', 'ARG', '91', '>10.4' ], [ 'I092R', 'ARG', '92', '>10.4' ], [ 'V099R', 'ARG', '99', '>10.4' ], [ 'N100R', 'ARG', '100', '>10.4' ], [ 'L103R', 'ARG', '103', '>10.4' ], [ 'V104R', 'ARG', '104', '>10.4' ], [ 'A109R', 'ARG', '109', '>10.4' ], [ 'N118R', 'ARG', '118', '>10.4' ], [ 'L125R', 'ARG', '125', '>10.4' ], [ 'A132R', 'ARG', '132', '>10.4' ] ]
    StatD = namedtuple('StatD', ['all_stat', 'glu_stat', 'asp_stat', 'his_stat', 'lys_stat', 'arg_stat'])
    extra_attrs = ('dsol', 'e_r','e_b','ele', 'ratio', 'dpk4', 'dpk20', 'ddpk')
        
    def __init__(self, session):
        temp = []
        for res in self.__class__.in_silico_res:
            if int(res[2]) < 50:
                res[2] = int(res[2]) - 6
            else:
                res[2] = int(res[2]) - 12
            #res[3] = float(res[3])
            temp.append(res)
        self.in_silico_res = temp
        self.dump = []
        self.ep = '4'
        self.session = session
        self.statds = ['all_dsol','low_dsol','mid_dsol','high_dsol']
        self.dsolranges = [(-999,999),(-999,3),(3,6),(6,999)]
        for statd in self.statds:
            all_stat = Stats('ALL Stats', self.__class__.extra_attrs)
            glu_stat = Stats('GLU Stats', self.__class__.extra_attrs)
            asp_stat = Stats('ASP Stats', self.__class__.extra_attrs)
            his_stat = Stats('HIS Stats', self.__class__.extra_attrs)
            lys_stat = Stats('LYS Stats', self.__class__.extra_attrs)
            arg_stat = Stats('ARG Stats', self.__class__.extra_attrs)
            self.__setattr__(statd, self.__class__.StatD(all_stat, glu_stat, asp_stat, his_stat, lys_stat, arg_stat))
    
    def Go(self):
        total = time.clock()
        now = time.clock()
        count = 0
        for prot in self.session.query(levels.ProtData).yield_per(10):
            rundict = {}
            fullrundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
                elif run.PACK!='f' and run.SCALE_VDW==.25:
                    fullrundict[run.EPSILON_PROT] = run
            if 4 in rundict and 20 in rundict and 4 in fullrundict and 20 in fullrundict:
                self.Kerchunk(fullrundict)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        print 'the whole thing took: %.3f' % (time.clock() - total)
        for item in self.dump:
            print item
            
    def Kerchunk(self, rundict):
        if 4 in rundict and 20 in rundict:
            for res4, res20 in zip(rundict[4].Ress, rundict[20].Ress):
                if res4==res20:
                    resep = eval('res'+self.ep)
                    res = resep.res.lower()
                    for dsolrange, statd in zip(self.dsolranges, self.statds):
                        if res in ['glu','asp','lys','arg'] and (dsolrange[0] < res4.dsol <= dsolrange[1]) and resep.dpk!=None and resep.seqadv=='': #\
                           #and (resep.Prot.pdb, resep.res, resep.resn) in zip(*zip(*self.in_silico_res)[0:3]):
                            for r in [res, 'all']:
                                if self.__getattribute__(statd).__getattribute__('%s_stat' % r).Gate([res4, res20]):
                                    if resep.dsol <= 0:
                                        dsol = 0.001
                                    else:
                                        dsol = resep.dsol
                                    if resep.dsol > 6 and res == 'asp':
                                        print resep
                                    self.__getattribute__(statd).__getattribute__('%s_stat' % r).dsol.append(dsol)
                                    self.__getattribute__(statd).__getattribute__('%s_stat' % r).e_r.append(resep.e_r)
                                    self.__getattribute__(statd).__getattribute__('%s_stat' % r).e_b.append(resep.e_b)
                                    self.__getattribute__(statd).__getattribute__('%s_stat' % r).ele.append(resep.e_r + resep.e_b)
                                    if res4.dsol != 0:
                                        self.__getattribute__(statd).__getattribute__('%s_stat' % r).ratio.append((resep.e_r + res4.e_b)/dsol)
                                    self.__getattribute__(statd).__getattribute__('%s_stat' % r).dpk4.append(res4.dpk)
                                    self.__getattribute__(statd).__getattribute__('%s_stat' % r).dpk20.append(res20.dpk)
                                    self.__getattribute__(statd).__getattribute__('%s_stat' % r).ddpk.append(res20.dpk - res4.dpk)
                else:
                    pass
    
    def Process(self):
        for statd in self.statds:
            print statd
            sd = self.__getattribute__(statd)
            for stat in sd:
                print stat

class EZPlotter(object):
    
    class Info(object):
        dsols = ('alldsol','lowdsol','middsol','highdsol')
        ress = ('ALL','ASP','GLU','ARG','LYS','HIS','TYR','SNASE')
        ResD = namedtuple('ResD', ('ALL','ASP','GLU','ARG','LYS','HIS','TYR','SNASE'))
        
        def __init__(self, **kwargs):
            D = namedtuple('D', ['stats']+kwargs.keys())
            for dsol in self.__class__.dsols:
                resargs = []
                for arg in self.__class__.ResD._fields:
                    resargs.append(None)
                self.__setattr__(dsol, self.__class__.ResD(*resargs))
                for res in self.__class__.ress:
                    dargs = [Stats('%s_%s' % (dsol, res))]
                    for val in kwargs.itervalues():
                        valargs = []
                        for arg in val._fields:
                            valargs.append([])
                        dargs.append(val(*valargs))
                    replacedict = {res:D(*dargs)}    
                    self.__setattr__(dsol, self.__getattribute__(dsol)._replace(**replacedict))
        
        def GenPlotInput(self, dsol, dname, x, y):
            plotinput = []
            for lis in (['ASP','GLU','ARG','LYS'], ['HIS'], ['TYR'],['SNASE']):
                for i in range(2):
                    plotinput.append([])
                for res in lis:
                    d = self.__getattribute__(dsol).__getattribute__(res).__getattribute__(dname)
                    plotinput[-2]+=d.__getattribute__(x)
                    plotinput[-1]+=d.__getattribute__(y)
            return plotinput
        
        def MakePlot(self, dname, outdir, outfile, plotter):
            plotinput = []
            settags = self.TagParse(outfile)
            for settag in settags:
                plotinput.append([])
                val,dsol,res = settag
                if res=='derk':
                    for subres in ['ASP','GLU','ARG','LYS']:
                        d = self.__getattribute__(dsol).__getattribute__(subres).__getattribute__(dname)
                        plotinput[-1]+=d.__getattribute__(val)
                else:
                    subress = []
                    for subres in ['ALL','ASP','GLU','ARG','LYS','HIS','TYR','SNASE']:
                        if subres.lower() in res:
                            subress.append(subres)
                    for subres in subress:
                        d = self.__getattribute__(dsol).__getattribute__(subres).__getattribute__(dname)
                        plotinput[-1]+=d.__getattribute__(val)
            plotter(*plotinput, outdir=outdir, outfile=outfile)
            
        @staticmethod
        def TagParse(outfile):
            rawtags = outfile.split('.')[-1].split('-')
            settags = []
            x = rawtags[0].split('_')[0]
            y = rawtags[0].split('_')[1]
            for dsol in rawtags[1].split('_'):
                for res in rawtags[2].split('_'):
                    settags.append((x, dsol, res))
                    settags.append((y, dsol, res))
            return settags

    houk_res = (('3RN3','ASP',14,1.8),('4LZT','ASP',87,2.1),('1PPF','ASP',27,2.2),('1XNB','ASP',11,2.5),('1BEO','ASP',21,2.5),('4LZT','ASP',18,2.7),('1XNB','ASP',106,2.7),('1PGA','ASP',22,2.9),('3RN3','ASP',121,3.1),('1A2P','ASP',75,3.1),('2RN2','ASP',94,3.2),('1PGA','ASP',47,3.4),('3RN3','ASP',53,3.7),('4LZT','ASP',52,3.7),('1PGA','ASP',36,3.8),('2TRX','ASP',20,3.8),('1DE3','ASP',59,4.1),('1DE3','ASP',57,4.3),('2RN2','ASP',10,6.1),('2TRX','ASP',26,8.1),('3RN3','GLU',2,2.6),('4LZT','GLU',7,2.9),('1PPF','GLU',19,3.2),('2RN2','GLU',57,3.2),('1A2P','GLU',60,3.2),('3RN3','GLU',111,3.5),('2RN2','GLU',129,3.6),('2RN2','GLU',61,3.9),('3RN3','GLU',9,4),('1BCA','GLU',26,4.1),('1PPF','GLU',10,4.1),('2RN2','GLU',119,4.1),('1PGA','GLU',27,4.5),('1PPF','GLU',43,4.8),('1DE3','GLU',96,5.1),('1ANS','GLU',20,5.4),('1RGA','GLU',28,5.9),('4LZT','GLU',35,6.2),('1HNG','GLU',41,6.7),('1XNB','GLU',172,6.7),('2SNM','LYS',66,6.4),('1L54','LYS',102,6.6),('1MUT','LYS',39,8.4),('1NFN','LYS',146,9.2),('1FEZ','LYS',53,9.3),('1GS9','LYS',146,9.4),('1LE2','LYS',143,9.4),('1NFN','LYS',143,9.5),('1NZP','LYS',312,9.5),('1GS9','LYS',143,9.9),('1LE2','LYS',146,9.9),('1PPF','LYS',34,10.1),('4LZT','LYS',33,10.4),('2BCA','LYS',41,10.8),('4LZT','LYS',96,10.8),('1PGA','LYS',28,10.9),('2BCA','LYS',16,10.9),('1PPF','LYS',55,11.1),('2BCA','LYS',7,11.2),('2BCA','LYS',55,12.4),('3EBX','HIS',6,2.8),('3SSI','HIS',43,3.2),('1STN','HIS',121,5.3),('4LZT','HIS',15,5.4),('1ERT','HIS',43,5.5),('1DE3','HIS',137,5.8),('3RN3','HIS',48,6.1),('3RN3','HIS',119,6.5),('3RN3','HIS',12,6),('1DE3','HIS',104,6.5),('1DE3','HIS',36,6.8),('2RN2','HIS',62,7),('2RN2','HIS',124,7.1),('1DE3','HIS',50,7.7),('1RGA','HIS',92,7.8),('1RGA','HIS',40,7.9),('2RN2','HIS',127,7.9),('1DG9','HIS',66,8.3),('2LZM','HIS',31,9.1),('1DG9','HIS',72,9.2))   
    in_silico_res = [ [ 'V023D', 'ASP', '23', '6.8' ], [ 'L025D', 'ASP', '25', '6.8' ], [ 'F034D', 'ASP', '34', '7.8' ], [ 'L036D', 'ASP', '36', '7.9' ], [ 'L037D', 'ASP', '37', '<4.0' ], [ 'L038D', 'ASP', '38', '6.8' ], [ 'V039D', 'ASP', '39', '8.1' ], [ 'T041D', 'ASP', '41', '<4.0' ], [ 'A058D', 'ASP', '58', '6.8' ], [ 'T062D', 'ASP', '62', '8.7' ], [ 'V066D', 'ASP', '66', '8.7' ], [ 'I072D', 'ASP', '72', '7.6' ], [ 'V074D', 'ASP', '74', '8.3' ], [ 'A090D', 'ASP', '90', '7.5' ], [ 'Y091D', 'ASP', '91', '7.2' ], [ 'I092D', 'ASP', '92', '8.1' ], [ 'V099D', 'ASP', '99', '8.5' ], [ 'N100D', 'ASP', '100', '6.9' ], [ 'L103D', 'ASP', '103', '8.7' ], [ 'V104D', 'ASP', '104', '9.7' ], [ 'A109D', 'ASP', '109', '7.5' ], [ 'N118D', 'ASP', '118', '7.0' ], [ 'L125D', 'ASP', '125', '7.6' ], [ 'A132D', 'ASP', '132', '7.0' ], [ 'V023E', 'GLU', '23', '7.1' ], [ 'L025E', 'GLU', '25', '7.5' ], [ 'F034E', 'GLU', '34', '7.3' ], [ 'L036E', 'GLU', '36', '8.7' ], [ 'L037E', 'GLU', '37', '5.2' ], [ 'L038E', 'GLU', '38', '7.0' ], [ 'V039E', 'GLU', '39', '8.2' ], [ 'T041E', 'GLU', '41', '6.5' ], [ 'A058E', 'GLU', '58', '7.7' ], [ 'T062E', 'GLU', '62', '7.7' ], [ 'V066E', 'GLU', '66', '8.5' ], [ 'I072E', 'GLU', '72', '7.3' ], [ 'V074E', 'GLU', '74', '7.8' ], [ 'A090E', 'GLU', '90', '6.4' ], [ 'Y091E', 'GLU', '91', '7.3' ], [ 'I092E', 'GLU', '92', '9.0' ], [ 'V099E', 'GLU', '99', '8.4' ], [ 'N100E', 'GLU', '100', '7.6' ], [ 'L103E', 'GLU', '103', '8.9' ], [ 'V104E', 'GLU', '104', '9.4' ], [ 'A109E', 'GLU', '109', '7.9' ], [ 'N118E', 'GLU', '118', '<4.5' ], [ 'L125E', 'GLU', '125', '9.1' ], [ 'A132E', 'GLU', '132', '7.0' ], [ 'V023K', 'LYS', '23', '7.3' ], [ 'L025K', 'LYS', '25', '6.3' ], [ 'F034K', 'LYS', '34', '7.1' ], [ 'L036K', 'LYS', '36', '7.2' ], [ 'L037K', 'LYS', '37', '>10.4' ], [ 'L038K', 'LYS', '38', '>10.4' ], [ 'V039K', 'LYS', '39', '9.0' ], [ 'T041K', 'LYS', '41', '9.3' ], [ 'A058K', 'LYS', '58', '10.4' ], [ 'T062K', 'LYS', '62', '8.1' ], [ 'V066K', 'LYS', '66', '5.6' ], [ 'I072K', 'LYS', '72', '8.6' ], [ 'V074K', 'LYS', '74', '7.4' ], [ 'A090K', 'LYS', '90', '8.6' ], [ 'Y091K', 'LYS', '91', '5.3' ], [ 'I092K', 'LYS', '92', '5.3' ], [ 'V099K', 'LYS', '99', '6.5' ], [ 'N100K', 'LYS', '100', '8.6' ], [ 'L103K', 'LYS', '103', '8.2' ], [ 'V104K', 'LYS', '104', '7.7' ], [ 'A109K', 'LYS', '109', '9.2' ], [ 'N118K', 'LYS', '118', '>10.4' ], [ 'L125K', 'LYS', '125', '6.2' ], [ 'A132K', 'LYS', '132', '>10.4' ], [ 'V023R', 'ARG', '23', '>10.4' ], [ 'L025R', 'ARG', '25', '>10.4' ], [ 'F034R', 'ARG', '34', '>10.4' ], [ 'L036R', 'ARG', '36', '>10.4' ], [ 'L037R', 'ARG', '37', '>10.4' ], [ 'L038R', 'ARG', '38', '>10.4' ], [ 'V039R', 'ARG', '39', '>10.4' ], [ 'T041R', 'ARG', '41', '>10.4' ], [ 'A058R', 'ARG', '58', '>10.4' ], [ 'T062R', 'ARG', '62', '>10.4' ], [ 'V066R', 'ARG', '66', '>10.4' ], [ 'I072R', 'ARG', '72', '>10.4' ], [ 'V074R', 'ARG', '74', '>10.4' ], [ 'A090R', 'ARG', '90', '>10.4' ], [ 'Y091R', 'ARG', '91', '>10.4' ], [ 'I092R', 'ARG', '92', '>10.4' ], [ 'V099R', 'ARG', '99', '>10.4' ], [ 'N100R', 'ARG', '100', '>10.4' ], [ 'L103R', 'ARG', '103', '>10.4' ], [ 'V104R', 'ARG', '104', '>10.4' ], [ 'A109R', 'ARG', '109', '>10.4' ], [ 'N118R', 'ARG', '118', '>10.4' ], [ 'L125R', 'ARG', '125', '>10.4' ], [ 'A132R', 'ARG', '132', '>10.4' ] ]
    
    DsolEleD = namedtuple('DsolEleD', ['dsol', 'ele', 'asite_dist'])
    RatioDDpkD = namedtuple('RatioD', ['ratio', 'ddpk'])
    DpkDDpkD = namedtuple('DpkDDpkD', ['dpk', 'ddpk'])
    
    Dpke4Dpke20D = namedtuple('Dpke4Dpke20D', ['dpke4', 'dpke20'])
    
    EleqElefD = namedtuple('EleqElefD', ['eleq', 'elef'])
    DsolqDsolfD = namedtuple('DsolqDsolfD', ['dsolq', 'dsolf'])
    DpkQDpkFD = namedtuple('DpkQDpkFD', ['dpkq', 'dpkf'])
    DDpkQDDpkFD = namedtuple('DDpkQDDpkFD', ['ddpkq', 'ddpkf'])
    
    NeuJiggleFEleAnom = namedtuple('NeuJiggleFEleAnom', ['neujigglef', 'eleanom'])
    NeuJiggleFDsolAnom = namedtuple('NeuJiggleFDsolAnom', ['neujigglef', 'dsolanom'])
    NeupKJiggleFEleAnom = namedtuple('NeupKJiggleFEleAnom', ['neupkjigglef', 'eleanom'])
    NeupKJiggleFDsolAnom = namedtuple('NeupKJiggleFDsolAnom', ['neupkjigglef', 'dsolanom'])
    IonJiggleFEleAnom = namedtuple('IonJiggleFEleAnom', ['ionjigglef', 'eleanom'])
    IonJiggleFDsolAnom = namedtuple('IonJiggleFDsolAnom', ['ionjigglef', 'dsolanom'])
    IonpKJiggleFEleAnom = namedtuple('IonpKJiggleFEleAnom', ['ionpkjigglef', 'eleanom'])
    IonpKJiggleFDsolAnom = namedtuple('IonpKJiggleFDsolAnom', ['ionpkjigglef', 'dsolanom'])
    
    def __init__(self, pdbsession=None, shsession=None, insilicosession=None, fullsession=None, csasession=None, consurfsession=None, outdir='/home/tel/benchmark/links'):
        temp = []
        for res in self.__class__.in_silico_res:
            if int(res[2]) < 50:
                res[2] = int(res[2]) - 6
            else:
                res[2] = int(res[2]) - 12
            #res[3] = float(res[3])
            temp.append(res)
        self.in_silico_res = temp
        self.dsoldict = {'alldsol':(-999,999),'lowdsol':(-999,3),'middsol':(3,6),'highdsol':(6,999)}
        
        self.pdbsession = pdbsession
        self.shsession = shsession
        self.insilicosession = insilicosession
        self.csasession = csasession
        self.consurfsession = consurfsession
        self.fullsession = fullsession
        self.outdir = outdir
        
        self.infokwargs = {}
        for key in ['DsolEleD', 'RatioDDpkD', 'DpkDDpkD', 'Dpke4Dpke20D', 'EleqElefD', 'DsolqDsolfD', 'DpkQDpkFD', 'DDpkQDDpkFD',
                    'NeuJiggleFEleAnom', 'NeuJiggleFDsolAnom', 'NeupKJiggleFEleAnom','NeupKJiggleFDsolAnom', 'IonJiggleFEleAnom', 'IonJiggleFDsolAnom','IonpKJiggleFEleAnom','IonpKJiggleFDsolAnom']:
            self.infokwargs[key.lower()] = self.__getattribute__(key)
        self.info = self.__class__.Info(**self.infokwargs)
    
    def Go(self):
        total = time.clock()
        now = time.clock()
        count = 0
        for prot in self.pdbsession.query(levels.ProtData).yield_per(10):
            rundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
            self.Kerchunk(rundict, snase=False)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        for prot in self.insilicosession.query(levels.ProtData).yield_per(10):
            rundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
            self.Kerchunk(rundict, snase=True)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        
        print 'graphing now...'
        self.info.MakePlot(dname='dsoleled', outdir=self.outdir, outfile='dsol_ele-alldsol-derk_his_tyr_snase-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='ratioddpkd', outdir=self.outdir, outfile='ratio_ddpk-highdsol-derk_his_tyr_snase-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='dpkddpkd', outdir=self.outdir, outfile='dpk_ddpk-alldsol-derk_his_tyr_snase-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='dpke4dpke20d', outdir=self.outdir, outfile='dpke4_dpke20-alldsol-derk_his_tyr_snase-multihexbin', plotter=MultiHexbinS)
        #SixPanelHisto(*(derk_ld+derk_md+derk_hd+snase_hd), outdir=self.outdir, outfile='dpka_ddpka-lowdsol_middsol_highdsol-derk-histos')
        #SixPanelHisto(*(his_ld+his_md+his_hd+snase_hd), outdir=self.outdir, outfile='dpka_ddpka-lowdsol_middsol_highdsol-his-histos')
        print 'the whole thing took: %.3f' % (time.clock() - total)
        
    def GoFQ(self):
        total = time.clock()
        now = time.clock()
        count = 0
        for prot in self.fullsession.query(levels.ProtData).yield_per(10):
            rundict = {}
            fullrundict = {}
            for run in prot.Runs:
                if run.PACK=='f':
                    rundict[run.EPSILON_PROT] = run
                elif run.PACK!='f' and run.SCALE_VDW==.25:
                    fullrundict[run.EPSILON_PROT] = run
            self.KerchunkFQ(rundict, fullrundict, snase=False)
            print '%d\t%s took: %.3f' % (count, prot.pdb, time.clock() - now)
            now = time.clock()
            count+=1
        
        print 'graphing now...'
        self.info.MakePlot(dname='eleqelefd', outdir=self.outdir, outfile='eleq_elef-alldsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='dsolqdsolfd', outdir=self.outdir, outfile='dsolq_dsolf-alldsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='dpkqdpkfd', outdir=self.outdir, outfile='dpkq_dpkf-alldsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='ddpkqddpkfd', outdir=self.outdir, outfile='ddpkq_ddpkf-alldsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)

        self.info.MakePlot(dname='neujigglefeleanom', outdir=self.outdir, outfile='neujigglef_eleanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='neujigglefdsolanom', outdir=self.outdir, outfile='neujigglef_dsolanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='neupkjigglefeleanom', outdir=self.outdir, outfile='neupkjigglef_eleanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='neupkjigglefdsolanom', outdir=self.outdir, outfile='neupkjigglef_dsolanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='ionjigglefeleanom', outdir=self.outdir, outfile='ionjigglef_eleanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='ionjigglefdsolanom', outdir=self.outdir, outfile='ionjigglef_dsolanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='ionpkjigglefeleanom', outdir=self.outdir, outfile='ionpkjigglef_eleanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        self.info.MakePlot(dname='ionpkjigglefdsolanom', outdir=self.outdir, outfile='ionpkjigglef_dsolanom-highdsol-derk_his_tyr-multihexbin', plotter=MultiHexbinS)
        print 'the whole thing took: %.3f' % (time.clock() - total)

    def Kerchunk(self, rundict, snase):
        if 4 in rundict and 20 in rundict:
            for res4, res20 in zip(rundict[4].Ress, rundict[20].Ress):
                if snase==True:
                    if (res4.Prot.pdb, res4.res, res4.resn) in zip(*zip(*self.in_silico_res)[0:3]):
                        pass
                    else:
                        continue
                if res4.dpk!=None and res4.seqadv=='':
                    if snase==True:
                        infolist = ['SNASE']
                    elif res4.res in ['ASP','GLU','ARG','LYS']:  
                        infolist = [res4.res,'ALL']
                    elif res4.res in ['HIS','TYR']:
                        infolist = [res4.res]
                    else:
                        infolist = []
                else:
                    continue
                for dsol, dsolrange in self.dsoldict.iteritems():
                    if (dsolrange[0] < res4.dsol <= dsolrange[1]):
                        for r in infolist:
                            info = self.info.__getattribute__(dsol).__getattribute__(r)
                            if info.stats.Gate([res4,res20]):
                                if res4.dsol<=0:
                                    dsol4 = .001
                                else:
                                    dsol4 = res4.dsol
                                info.dsoleled.dsol.append(dsol4)
                                info.dsoleled.ele.append(res4.e_b + res4.e_r)
                                info.dpkddpkd.dpk.append(res4.dpk)
                                info.dpkddpkd.ddpk.append(res20.dpk - res4.dpk)
                                info.dpke4dpke20d.dpke4.append(res4.dpk)
                                info.dpke4dpke20d.dpke20.append(res20.dpk)
                                if dsol4 > .5:
                                    info.ratioddpkd.ratio.append((res4.e_b + res4.e_r)/dsol4)
                                    info.ratioddpkd.ddpk.append(res20.dpk - res4.dpk)
    
    def KerchunkFQ(self, rundict, fullrundict, snase):
        '''
        kerchunker for full vs quick analysis
        '''
        if 4 in rundict and 20 in rundict and 4 in fullrundict and 20 in fullrundict:
            for qres4, qres20, fres4, fres20 in zip(rundict[4].Ress, rundict[20].Ress, fullrundict[4].Ress, fullrundict[20].Ress):
                if snase==True:
                    if (qres4.Prot.pdb, qres4.res, qres4.resn) in zip(*zip(*self.in_silico_res)[0:3]):
                        pass
                    else:
                        continue
                if qres4==qres20==fres4==fres20 and qres4.dpk!=None and qres4.seqadv=='':
                    if snase==True:
                        infolist = ['SNASE']
                    elif qres4.res in ['ASP','GLU','ARG','LYS']:  
                        infolist = [qres4.res,'ALL']
                    elif qres4.res in ['HIS','TYR']:
                        infolist = [qres4.res]
                    else:
                        infolist = []
                else:
                    continue
                for dsol, dsolrange in self.dsoldict.iteritems():
                    if (dsolrange[0] < qres4.dsol <= dsolrange[1]):
                        for r in infolist:
                            info = self.info.__getattribute__(dsol).__getattribute__(r)
                            if info.stats.Gate([qres4,qres20,fres4,fres20]):
                                info.eleqelefd.eleq.append(qres4.e_b + qres4.e_r)
                                info.eleqelefd.elef.append(fres4.e_b + fres4.e_r)
                                info.dsolqdsolfd.dsolq.append(qres4.dsol)
                                info.dsolqdsolfd.dsolf.append(fres4.dsol)
                                info.dpkqdpkfd.dpkq.append(qres4.dpk)
                                info.dpkqdpkfd.dpkf.append(fres4.dpk)
                                info.ddpkqddpkfd.ddpkq.append(qres20.dpk - qres4.dpk)
                                info.ddpkqddpkfd.ddpkf.append(fres20.dpk - fres4.dpk)
                                
                                neujigglef = fres4.GetJiggle()
                                neupkjigglef = fres4.GetJiggle(atpk=True)
                                ionjigglef = fres4.GetJiggle(ionized=True)
                                ionpkjigglef = fres4.GetJiggle(ionized=True, atpk=True)
                                
                                if ionjigglef < 1:
                                    info.neujigglefeleanom.neujigglef.append(neujigglef)
                                    info.neujigglefeleanom.eleanom.append((fres4.e_b + fres4.e_r) - (qres4.e_b + qres4.e_r))
                                    info.neujigglefdsolanom.neujigglef.append(neujigglef)
                                    info.neujigglefdsolanom.dsolanom.append(fres4.dsol - qres4.dsol)
                                
                                if ionpkjigglef < 1:
                                    info.neupkjigglefeleanom.neupkjigglef.append(neupkjigglef)
                                    info.neupkjigglefeleanom.eleanom.append((fres4.e_b + fres4.e_r) - (qres4.e_b + qres4.e_r))
                                    info.neupkjigglefdsolanom.neupkjigglef.append(neupkjigglef)
                                    info.neupkjigglefdsolanom.dsolanom.append(fres4.dsol - qres4.dsol)
                                
                                info.ionjigglefeleanom.ionjigglef.append(ionjigglef)
                                info.ionjigglefeleanom.eleanom.append((fres4.e_b + fres4.e_r) - (qres4.e_b + qres4.e_r))
                                info.ionjigglefdsolanom.ionjigglef.append(ionjigglef)
                                info.ionjigglefdsolanom.dsolanom.append(fres4.dsol - qres4.dsol)
                                
                                info.ionpkjigglefeleanom.ionpkjigglef.append(ionpkjigglef)
                                info.ionpkjigglefeleanom.eleanom.append((fres4.e_b + fres4.e_r) - (qres4.e_b + qres4.e_r))
                                info.ionpkjigglefdsolanom.ionpkjigglef.append(ionpkjigglef)
                                info.ionpkjigglefdsolanom.dsolanom.append(fres4.dsol - qres4.dsol)
        
    def Process(self):
        for stat in self.statd:
            print stat
            
if __name__=='__main__':
    dbpath = sys.argv[1]
    
    dbdict = {}
    metadict = {}
    sessiondict = {}
    #dbl = ['in_silico']
    #dbl = ['dunbrack','sh','actprot', 'in_silico', 'full']
    dbl = ['full']
    #dbl = ['foo']
    
    for dbn in (dbn for dbn in dbl if os.path.isfile(os.path.join(dbpath, dbn+'.db'))):
        dbdict[dbn] = create_engine('sqlite:///' + os.path.join(dbpath, dbn+'.db'), poolclass=NullPool)
        event.listen(dbdict[dbn], 'connect', _fk_pragma_on_connect)
        metadict[dbn] = MetaData(dbdict[dbn])
        
        for level in levels.hierarchy:
            level.Table(meta=metadict[dbn], db=dbdict[dbn], reloadt=True)
        sessiondict[dbn] = LoadSession(dbdict[dbn])
        
    if os.path.isfile(os.path.join(dbpath,'csa.db')):
        csadb = create_engine('sqlite:///' + os.path.join(dbpath,'csa.db'), poolclass=NullPool)
        event.listen(csadb, 'connect', _fk_pragma_on_connect)
        csameta= MetaData(csadb)
        levels.CSAData.Table(meta=csameta, db=csadb, reloadt=True)
        csasession = LoadSession(csadb)
    
    graphs = EZPlotter(fullsession=sessiondict['full'], csasession=csasession, outdir=dbpath)
    #sessiondict['dunbrack'], sessiondict['sh'], sessiondict['in_silico'],
    graphs.GoFQ()
    #graphs.GoSHIndividual()
    graphs.Process()
    
#    averages = Averages(session=sessiondict['full'])
#    averages.Go()
#    averages.Process()